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Abstract 

Using the specific model of a system of like charged ions confined between two 
planar like charged surfaces, we compare the predictions for the energy and 
density profile of four simulation methods available to treat the long range 
Coulomb interactions in systems periodic in two directions but bound in the 
third one. Monte Carlo simulations demonstrate unambiguously complete 
agreement between the results obtained with these methods where the poten- 
tial between charges is solution of Poisson's equation in the simulation cell 
with adequate boundary conditions. The practical advantages of the different 
methods is assessed, imulations, Monte Carlo, Coulomb interaction, colloids 
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I. INTRODUCTION 

The standard way to treat long range interactions (Coulomb, dipolar, Yukawa near the 
Coulomb limit) in simulations of bulk systems is to replicate the basic simulation cell peri- 
odically in all three directions of space and apply an Ewald summation technique (EW3D) 
imj^. However, in many situations of electrochemical, biological or technological interest, 
as, for instance, electrolyte solutions between charged surfaces, charged lipid bilayers in wa- 
ter, suspensions of colloids between glass plates, clays magnetic thin films 0, Wigner 
crystals 0] etc., the system is finite in one direction thus necessitating a modification of the 
convential Ewald method. For systems periodic in two dimensions but bounded along the 
third one an Ewald summation method (EW2D) for electrostatic interactions has first been 
given by Parry p],^ and was later rederived by various authors [p!0|-p!^. Unfortunately, the 
practical use of the EW2D sum is hampered by the occurrence, in the reciprocal space term, 
of a double sum over different particles which, due to the complicated way the bounded 
coordinates enter the expression, cannot be reduced, as for the EW3D case, to a sum of 
order N, a circumstance which renders the method computationally expensive. Not sur- 
prisingly, only few calculations using EW2D have been reported to date mainly to test the 



validity of more approximate approaches • These calculations involve tabulation and 



inteerpolation of the potentials on a three-dimensional grid. 

Various "time saving" schemes have been proposed to bypass the computational burden 
of the EW2D method. The purpose of this paper is to compare some of these methods 
for the specific case of a system of charged particles confined between two charged planar 
surfaces separated by a distance of the order of several particle diameters. The methods 
that have been chosen are those which satisfy exactly the laws of electrostatics (perfect 
screening. Green's theorem etc..) and do not present numerical instabilities difficult to 
control. Thus the charged sheet method proposed by Torrie and Valleau [|T3)|I1] ^^^1 its 
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modification by Boda et al [|T9[ have been discarded. Indeed the distribution of the point 
charges of the ions outside the simulation cell, located in the image cells of the latter, is 
represented approximately by a set of uniformely charged planar sheets. An other method, 
proposed by Lekner [0,^ is based on rewriting the sum of forces acting on an ion by a 
second ion and its periodic images as an absolutely converging infinite sum. The energy is 
obtained by integration of the force. A shortcoming of the method is the need, to estimate 
the sum with a given precision, to retain a number of terms depending strongly on the 
relative position of pairs of particles. Because of this technical reason, emphasized in more 
detail in ref [^, the method was not included in this investigation. 

The four methods we have considered, the first three of which have already appeared 
in the literature, are EW3D |]r5| , p!6| , p3| , p4| , a method developed by Hautman and Klein ||2^ 



(HK), the method of hyperspheres (HSG) and that of concentric spheres (CS) described 
below. 

A way to treat the long range Coulomb interactions is to place the slab at the centre 
of a parallelepipedic simulation box which has dimension perpendicular to the slab much 
larger than the width of the latter and apply the EW3D method. In this way one effectively 
simulates an infinite number of parallel slabs. Provided the region of empty space separating 
the slabs is sufficiently large one expects the influence of the periodic images on the behaviour 
of the system to become negligible. 

If the distance between the surfaces confining the particles is small compared to the 
extension of the system along the directions parallel to the surfaces (narrow slab), the 
lateral distance s between pairs of particles will be much larger than the distance z normal 
to the slab surfaces and it will be appropriate to expand the Coulomb pair interaction in 
powers of z/s and apply an Ewald sum to the in-plane component of the interaction. Such 
an approach was followed by Hautman and Klein (HK) ||25|| . 

An attractive way to investigate the behaviour of Coulomb particles between charged 
surfaces which circumvents the cumbersome Ewald sums is to use a closed space, e. g. a 



hypersphere in four- dimensional Euclidean space [Pq -|2q|. Not only is the electrostatic po- 



tential, solution of Poisson's equation on the hypersphere obtainable in simple closed form 
p7| , thus avoiding approximation of the interaction potential (as for instance in EW3D due 
to truncation of the direct and reciprocal space sums), also the number of operations neces- 
sary to calculate the distances between particles is reduced with respect to Euclidean space 
with concomitant speed up of the simulation procedure. This geometry has been applied 
previously to the study of the attraction between two like charged surfaces neutralized by 
solvated counterions in an endeavour to comprehend the stability of charged lamellar mate- 
rials such as clays and cement and to the study of the effective interaction of charged 



colloidal particles confined between like charged plates ||30|1 . 

A system of charged particles occupying the region between the outer and inner surfaces 
of two concentric spheres is also a suitable arrangement to describe, in the limit of sufficiently 
large radii of the spheres, the system of charged particles confined between two planar 
surfaces. This geometry has the advantage that the interaction between charges and between 
charges and surfaces is the usual Coulomb potential; however it may require use of a large 
number of particles to render the effect of curvature of the confining surfaces small. 

The purpose of this paper is to check, by means of Monte Carlo simulation, the agreement 
between the results of the four aforementioned methods for the energy and density profile of 
the system described above consisting of N like-charged ions between parallel like-charged 
surfaces separated by a distance h. The ions, modelled as hard spheres of diameter d 
bear a charge q at their centre and each surface, of area S, a uniform charge density a. 
Such a system can be viewed as a crude model for lamellar liquid crystals formed by ionic 
amphiphiles |^ or charged lamellar materials like clays or cement [|2^. The convergence 



of the results to their thermodynamic limit is estimated by performing simulations with an 
increasing number of particles keeping the distance between surfaces fixed. The speed of 
this convergence is obviously an important criterion for appraisal of the relative practical 
interest of the four methods. 

Expressions for the energy calculated with the different methods will be given in Sect. 
II together with details of the Monte Carlo simulations. Results for the energy and density 



profiles obtained with the four methods are compared in Sect. Ill and conclusions drawn in 
Sect. IV. 



II. ENERGY EXPRESSIONS 



A. EW3D 



In this method a square slab of ions of thickness h is placed at the centre of a simulation 
box having dimension normal to the surfaces much larger then the lateral dimension L 
and the system is extended periodically in the three directions of space. The slab surfaces 
are located z — ±h/2 perpendicular to the 2;-axis. The closest approaches of the ions to 



To evaluate the total energy of the system which is the sum of the contributions from the 
ion-ion Ua, ion-surface Uig and surface-surface Ugs interactions it will be convenient to asso- 
ciate to the ions a uniform neutralizing background (filling the whole simulation cell) of den- 
sity —Nq/V and to the charged walls a uniform background of density —2aS/V — —2a/Lz 
where V — L^L'^ — L^S is the volume of the simulation cell. Because of electroneutrahty of 
the system Nq + 2a S = 0, the backgrounds cancel out exactly and will not contribute to 
the total energy. 

The contribution of the ions and their background to the total energy is given by 



In these equations Vij — Vj — Vi where rj and Vj are the positions of particles i and j. 



the impenetrable surfaces are therefore z = ±^ 



(h-d). 




(1) 



The term —\Cy, is the self-energy given by 
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2« TT I 

+ ^ + 7^ (2) 



\/vr . 

In Eq. (|I|) 1/ is a vector of components {Lnx,Lny,Lznz) {nx,ny,nz integers) and the prime in 
the sum over v means that the terms i = j must be omitted when u = 0. The wave- vectors 
G which enter the reciprocal space contributions to the energy have components 27m;j./L, 
2TTny/L and 27mz/Lz 

In our calculations the sum on lattice vectors extends over all G subject to \nx\ < 6, 
\ny\ < 6, \nz\ < 12 and |n| < rimax = 12. The a parameter which governs the rate of 
convergence of the real- and reciprocal-space contributions was taken sufficiently large so 
that only the terms with u = had to be retained in Eqs. (|l|) and 

The electrostatic energy between the ions and the two charged walls including their 
compensating backgrounds is given by 

= E qVs{z,) - ^ / drVsiz) (3) 
i=i ^ "'^ 

where Vs{z) is the electric potential associated with the electric field Eg generated by 
the two charged surfaces and theit backgrounds. Due to the symmetry of the simula- 
tion cell, this electric field is normal to the charged surfaces and depends only on the z- 
coordinate. Furthermore, as a consequence of the periodic boundary conditions, it must 
satisfy Es{—Lz/2) = Es{Lz/2) = 0. Applying Gauss's theorem (cf. Fig. 1) one finds 

TT(J 

z - ina -Lz/2 <z< -h/2 



Es{z) 



z 



L 

-^z -h/2 <z <h/2 (4) 

Lz 

'^^-z + Ana h/2 <z < Lz/2 



Lz 

and, by integration, for the potential, choosing Vs{0) = (any additional constant would 
leave Uis unchanged). 
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VJz) 



+ 4:7raz + 2nah -LJ2 < z < -h/2 



Aira 



Ana 



z — Airaz + 27iah 



-h/2 <z <h/2 
h/2 <z < LJ2 



(5) 



In the periodic system the potential created by the charged surfaces is thus parabohc. The 
potential Vs{z) can also be obtained by a direct integration of the Ewald potential over 
the two charged planes resulting in a one-dimensional Fourier series which can be summed 
explicitly to yield Eq. Integration of K over the volume of the simulation box leads to 



N 

E 



Z; — 



Nqna 



(3/1^ - 6Lzh + 2L 



(6) 



Finally, the interaction between the two charged surfaces (including their background) is 

2nScx^ 



-[h' + {Lz-hY]. 



(7) 



B. Hautmann-Klein method. 

In this method the slab, having the characteristics described above, has periodic bound- 
ary conditions only in the x- and y- directions. The origin of coordinates being at the centre 
of the parallelepipedic cell, the energy of the system is given by 

^2 N 1 Af 2 



V E E '1 ; — f + E E / d^adyc. E 1 : — r 

2 ,.-7^1 „ I r,-,; + u \ ~i ^Js „ I r,; - r„ + 1/ I 



U 

2 



i,j=l 1/ I ' «J I " I i=l a=l " ■5' U 



Iff 1 

+ 2^' E / dx'Jy'^ dxpdyp E i _ ^ , ^ 



■.13= 

where Vij = Vj — r^, ri,rj coordinates of the particles i and j, and the vectors Va and 
are the position vectors of points on the surfaces Sa (a = 1, 2) having constant z coordinate 
z = h/2 if a or /? = 1 and z = —h/2 if a or /3 = 2. The sum over u runs over all vectors of 
components {n^L, riyL) (n^., Uy integers) perpendicular to the z-direction. The three sums 
in the expression for U correspond to the interactions between the particles, the interaction 
between the particles and the surfaces Si and 5*2 and to that between the two surfaces 5*1 
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and 5*2, respectively. To each of these sums, due to the infinite number of rephcated cells, 
is associated a divergent contribution. However, if electroneutrality is taken into account 
these divergent terms will cancel. 

The evaluation of the ion-ion interaction starts with the identity |]25[ 



- = (- - E (^n-^) + E ('n-^ (9) 
' ' n=0 n=0 

with 

_ (-l)"(2n)! 
22"(n!)2 ■ 

The added and subtracted term represents the first M + 1 terms in the binomial expansion 
of 1/r in powers of z/s where s is the component of r in the plane of the surface. By 
introducing a convergence function /i„(s; a) for each term l/s^""*"^ the energy can be divided 
in a short range part 

-|A''A'' -| M f \ 



2 i=i j=i D fij,!^ 



V ' n=0 ^ijM 



and a long range part 

T N N M h { s ■ 

ul = Ij:j: E «^4"( E (n) 

^ j=l j=l n=0 V ^ij,U 

which can be evaluated in reciprocal space. In these equations Sij^iy = \sj — Si + v\. With 
the choice 

ho{s; a) = erf(as) (12) 

and 

hn{s;a) _ 1 y,2n( ho{s;a) \ 

a„(2n)! ^ s ^ ^ ^ 

as proposed by Hautman and Klein the short and long ranged parts of the electrostatic 



energy take the form [25 



^ i=l jyi ' tj ra=l v^'v- 



where only the u = term has been retained (see below) and 
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E E m E 

i=l j=l n=0 



(2n) 



2n 



(15) 



E^.-^(E^O 



with G a two-dimensional vector of components 27r(n^./L, Uy/ L). In our simulations we kept 
terms up to M = 3 which allowed the short range part f/s to be limited to its = term 
even for relatively large values of h. Expressions for /i„(s; a) up to M=3 are 



hi{s]a) = erf(Q;s) —e 



"'^Vl + 2aV) 



h2{s] a) = erf(as) — 



2as 



e-"'^'(l + -aV 



2 
-( 
3 



^44 

-as 



(16) 



hsis; a) 



erf(as) 



2as 



3 15 25 



225 225 ^ 

The attractive feature of the method is that, by writing the zfj' explicitly as polynomials 
in Zi and zj, U-^ becomes a sum of terms each of which involvs products of two functions 
having general form 



(17) 



In the evaluation of U^^ sum on the pairs of particles is replaced by sums on particles which, 
obviously, makes the computation faster. The contributions to the energy of the interactions 
between the ions and the surfaces and between the surfaces are easily obtained using the 



method of de Leeuw and Perram ||Tl]] and the identity 



r /'+00 /'+00 

dxa,dya,^f{\ra-r + u\)= dx^ dy^ f{\ Va - r \). 



(18) 
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The divergent terms of these two contributions having been ehminated through use of the 
electroneutrahty condition, the ion-surface and surface-surface interactions contribute to the 
energy f/ by a constant term equal to 27ra'^V. The expression for the total energy is therefore 
given by 

f/ = f75 + f/^ + 2na'V. (19) 



C. Hyperspherical geometry 

In this method the Monte Carlo simulations are performed on a hypersphere 5*3 in four- 
dimensional Euclidian space. Two surfaces of angular colatitudes 9n and 9s = n — Ojsi 
separated by a distance h = R{tt — 26 n) are located symmetrically on opposite sides of the 
equator (see Fig. 3 of ref. N neutralizing ions of charge q are confined between the two 

surfaces which bear each a charge density a. For given h and a the aperture 6' at is obtained 
by solving the equation 

hsme^ ( qN^' 

n-2ej,~[ 8na) ' ^^^^ 
Similar to the EW3D case described above, neutralizing backgrounds are associated to 
both the ions and the charged surfaces. As a detailed derivation of the different contributions 
to the potentiel energy, ion-ion Uu, ion-surface Uis and surface-surface Uss can be found in 
ref. only the final expressions are given here 



2 Af A'' 

= ^ E Ei(^ - - 0.5] 

H ^ [(tt - 2a)cota + 0.5 - 7r/sina] (21) 

Uis = ^[(vr - 2^i)cot^, - 1] + ^(^jvcot^Tv - 1) (22) 
Uss = \[l + {'r^-^9N)cot9j,]. (23) 

Here the angle a is equal to the ion radius d/2 divided by the hypersphere radius i?, Oij is 
the angular separation between particles i and j and = o^S {S area of each surface). 
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D. Concentric spheres 



In the last method we have explored, the ions occupied the region of volume V between 
two concentric spheres of radii ri and = h + ri. The external surface of the inner sphere, 
of area Si = 4:7Trf, and the internal surface of the outer sphere, of area Sm = 47rr^, bear a 
charge density a. The two surfaces being impenetrable the distances of closest approach of 
an ion to the surfaces are ri + d/2 and Vm — d/2. The electroneutrality conditions reads 

a{Si + Sm) + Nq = 0. (24) 

The total energy of the system is readily obtained as 

where the three first terms in the r.h.s. represent the ion-ion energy, the energy of the 
ions with the charged surface Si and with the surface 5*^ which creates in the volume V 
a constant potential. The three last terms correspond to the energies associated with the 
surface charges of Si and Sm- 



III. RESULTS 

In our comparisons we did not consider realistic values of q and a as for instance envisaged 
in ref. ||27|j29|l . Our aim being methodological, we chose values of q and cr allowing for a 



compelling test of the efficiency of the four methods to predict reliably the properties of the 
confined charged system. In the following charges q and a will be expressed by using reduced 
units of charge and distance and d, respectively. In these units the surface charge 

a has been fixed to -1 in all our simulations and the value of q varied from 0.5 to 5. The 
distance between the surfaces limiting the system has been taken to be /i = 5; however, 
in order to test the limit of validity of the HK method for large h a few simulations were 
performed with h = 8 and h = 12. For given values of h and a electroneutrality entails that 
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the density of the system increases for decreasing values of q if the volume of the simulation 
cell remains unchanged. 

The simulations were carried out in the canonical ensemble. Of the order of 10^ MC steps 
per particle were performed to calculate the energy and the density profiles when < 3000. 
In the CS geometry, for 15000, 10^ trial configurations per particle were generated 

to calculate these quantities. The statistical error on the energies is of the order of ±0.02 
in units of kT. For the two planar geometries, EW3D and HK, the density profiles are 
estimated with an error of ~ 0.5%; the statistical error, except for systematic error due to 
curvature effects, is of the same order for the HSG method. In the case of the CS geometry, 
for 15000, the error on the energies is of the order of ±0.03 and on the density profiles 
of the order of 2 — 3%. 

In the calculations using the EW3D method, the simulation volume has a square section 
of side L and an elongation along the 2;-axis of varying between 60 and 90. For h = 5 and 
L between 15 and 30, the influence of the value of on the simulation results turned out to 
be always negligible. Also, within this geometry, if the system confined between the surfaces 
has a net dipole moment, a correction to the energy proportional to the square of the dipole 
moment normal to the surfaces can be taken into account to remove the interaction between 



net dipoles in the periodically repeated slabs ||16| . As the dipole moment in our system was 
small such a correction was not considered. 

A systematic comparison of the energies for 0.75 < g < 5 is made in Table 1. It shows 
without ambiguity the convergence of the results of the four methods. It appears that the 
thermodynamic limit for U is obtained for g > 1 as soon as the number of particles is of the 
order of 500 for EW3D and HK and of the order of 2000-3000 for the HSG and CS method. 
However, at g ^ 0.75, i. e. at densities larger than 0.4 at fixed h and a, at least 10000 
particles are necessary for the energy calculated with CS to agree within 1% with the other 
methods. This result is obviously related to the important curvature effects expected for 
radii r; and < 10. Figure |^ shows, for the four methods, the variation of U when the ionic 
charge q varies from 0.5 to 5. A maximum is observed for q ^ 0.8. This maximum occurs as a 
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consequence of the lowering of q and saturation of the particles in contact with the confining 
surfaces due to steric effects. Indeed layers of particles parallel to the surfaces form at high 
density as evidenced by Fig. |^ which presents the density profiles for the different values of 
q. For large g, layers of highly localized particles are in contact with the surfaces. When 
the density increases (i. e. q decreases) layers separated by a distance of order d appear 
in the volume. The profiles shown in Fig. |^ are obtained by the method of HK and are 
indistinguishable from those given by EW3D. 

Figures § and ^ allow to compare the density profiles obtained with the methods EW3D 
and HK with those evaluated with the methods of hypersphere and concentric spheres. For 
q = 2 excellent agreement is found for the four methods taking into account the statistical 
errors , especially those for the CS method. At q = 0.75 differences are manifest and 
demonstrate the persistence of curvature effects; they are small in the hypersphere geometry 
but remain important in the CS geometry in spite of the large value N = 14036. This result 
at g = 0.75 is in agreement with the slow convergence, in the latter geometry, of the energy 
to its thermodynamic limit when g < 1. 

Extrapolations towards z = ±2 give the values pc of the profiles when the particles are 
in contact with the surfaces. It has been shown that for planar surfaces separated by a 
distance much larger than the diameter d of the particles the pressure is given by |32| , |33|| 

P/kT = - El/MkT (26) 

where Ec-, the value of the electric field at the surface, equals Aira for planar surfaces. For 
the confined systems considered in this work this expression is not strictly valid anymore. 
Nonetheless, considering it as an approximation and using the values of pc given in Table 2, 
an estimate of P can be obtained. The value of P is small and positive for g > 1 and rises 
to a value of 5 at g = 0.5 where hard core interactions are important. The nearly zero value 
of the pressure for g > 1 is clearly in accord with the absence of particles in the central 
part of the simulation cell (cf. Fig. ^ and complete screening of the charged surfaces by the 
particles in their vicinity implying an almost vanishing electric field in this same region of 
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simulation volume. 

Table 2 and Fig. ^ comparing simulation results for the energy and density profiles for 
h = 8 and 12 obtained with the EW3D and HK methods, show that the latter method 
remains accurate even for /i/L ~ 0.5. However, reliable results are obtained only if, in Eqs. 
H| and [T^ for the energy, M is chosen to be 3. 



IV. CONCLUSION 

The convergence of the results obtained with the four simulation methods considered 
is clearly the main conclusion of this study. It demonstrates the strict equivalence of the 
various possible routes to take into account the long range of the Coulomb interactions : 
use of periodic boundary conditions, a closed system without boundary or simply a large 
system with boundaries. The equivalence is made realized only by the use of potentials 
preserving the laws of electrostatics and thus solution of Poisson's equation associated with 
the simulation cell and its boundary conditions. One can remark that the system under 
investigation is particularly challenging to establish this equivalence because particles all 
bear a charge of the same sign and screening effects therefore result only from interactions 
between the particles and the surfaces. 

Depending on the type of system which is simulated, the advantage of using one method 
or the other may differ. This work clearly shows that the EW3D and HK methods apply 
more favourably to the study of systems confined by planar surfaces, the thermodynamic 
limit being reached for ~ 500. However, with regard to computing time, the method of 
hypersphere is the most efficient despite the necessity to use systems with larger values of 

to make curvature effects negligible. Obviously, the CS method is unfavourable to study 
planar interfaces as curvature effects get small only for N > 10000 — 20000. It has been 
used here mainly to establish without ambiguity the equivalence between the Ewald and 
hypersphere "Coulomb" potentials and the usual Coulomb potential. 

The calculations carried out with the CS method nevertheless show that the ionic density 
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profiles are affected in the vicinity of a charged interface by the curvature of the latter. 
Such modifications of the density profiles near planar charged interfaces may be present in 
solutions of inverted micelles or suspensions of charged colloids of small diameter. 

The determination of the relative efficiency of the simulation methods for the calculation 
of the properties of confined systems will be complete only when pressure, free energy and 
surface tension will have been evaluated. Calculation of these quantities is in progress. 
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FIGURES 

FIG. 1. Simulation setup for EW3D. The slab of width h is placed at the center of a simulation 

cell having dimension in the direction normal to the slab larger than the lateral dimensions 
L. Periodic boundary conditions are applied in all three directions. To evaluate the electric field 
Es{z) created by the charged surfaces and their backgrounds, Gauss's theorem is applied to the 
dotted volume. 

FIG. 2. Variation with ionic charge q of the total energy of a system of like ions confined 
between planar surfaces bearing a charge density a = — 1 and separated by a distance h = 5, for 
the four different simulation methods considered. 

FIG. 3. Variation with ionic charge q of the density profile p{z) of a system of like ions confined 
between planar surfaces bearing a charge density a = —1 and separated by a distance h = 5. The 
results shown correspond to the method of HK and are indistinguishable from those obtained for 
EW3D. 

FIG. 4. Ion density profile for a slab of width h = 5, surface charge density a = —1 and ionic 
charge q = 2. Comparison between all four simulation methods. The profile is shown for z/d > 1. 

FIG. 5. Ion density profile for a slab of width h = 5, surface charge density a = —1 and 
ionic charge q = 0.75. Comparison between all four simulation methods. The profile is shown for 
z/d > -0.5. 

FIG. 6. Comparison of the density profiles obtained with the HK and EW3D methods for (from 
top to bottom) h = 5, a = —1; h = 8, a = —1 and h = 12, a = —2. All results are for q = \. For 
clarity the profiles have been shifted by 0.5 with respect to each other. 
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TABLES 

TABLE I. Variation with ionic charge q of the total energy of a system of like ions confined 

between planar surfaces bearing a charge density a = —1 and separated by a distance h = 5, 
for the four different simulation methods considered. The number N of ions is given in brackets, 
P = l/kT, T temperature. 





/3Uew3d/N 


PUhk/N 


PUhsg/N 


PUcs/N 


5.0 


-4.30 (320) 
-4.30 (980) 


-4.31 (1024) 


-4.26 (1024) 
-4.29 (3000) 


-4.16 (1024) 
-4.26 (3544) 
-4.28 (14036) 


4.0 


-1.30 (320) 
-1.28 (500) 


-1.28 (1024) 


-1.26 (1024) 
-1.27 (3000) 




3.0 


0.930 (500) 


0.943 (1024) 


0.948 (1024) 

0.94G (3000) 




2.0 


2.32 (500) 
2.35 (980) 


2.35 (1024) 


2.36 (1024) 
2.35 (3000) 


2.33 (1024) 
2.35 (3544) 

2.3G (11036) 


1.0 


3.32 (980) 


3.32 (1024) 


3.23 (1024) 
3.29 (3000) 


3.51 (1024) 
3.36 (3544) 

3.30 (14036) 


0.75 


3.13 (500) 
3.17 (980) 
3.16 (1620) 


3.15 (1024) 


3.02 (1024) 
3.10 (3000) 


3.57 (1024) 
3.26 (3544) 
3.19 (14036) 
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TABLE II. Comparison of energy U and contact value pc of the density profile for the two 
methods HK and EW3D. The surface charge density is cr = —1 except for the case h = 12 where 
a = —2, /? = 1/kT, T temperature. The average density is p = N/V = —2a/qh. 





Ewald 3d 


Hautman-Klein 


tl 


q 


P 


IDTT 1 'AT 

pU /IS 


Pc 


QJJ 1 AT 

pU /T\ 


Pc 


5.0 


5.0 


0.08 


-4.30 


6.30 


-4.31 


6.38 




4.0 


0.1 


-1.28 


6.30 


-1.28 


6.32 




3.0 


2/15 


0.93 


6.30 


0.94 


6.31 




2.5 


0.16 


1.70 


6.30 


1.71 


6.29 




2.0 


0.2 


2.35 


6.30 


2.35 


6.29 




1.5 


4/15 


2.97 


6.35 


2.98 


6.33 




i.U 




o.oZ 


D.Oo 








0.75 


8/15 


3.17 


7.22 


3.15 


7.20 




0.625 


0.64 


2.90 


8.36 


2.90 


8.38 




0.5 


0.8 


2.49 


11.6 


2.49 


11.6 


8.0 


1.0 


0.25 


3.46 


6.40 


3.53 


6.39 


12.0" 


1.0 


1/3 


9.07 


27.4 


8.92 


29.3 



V = -2 
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/K 





1 






L 


. J 








h 

2 



h 
2 



